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Abstract 

The aeroelastic characteristics of a compressor stage were analyzed using a computational fluid 
dynamic (CFD) solver that uses the harmonic balance method to solve the governing equations. The three 
dimensional solver models the unsteady flow field due to blade vibration using the Reynolds-Averaged 
Navier-Stokes equations. The formulation enables the study of the effect of blade row interaction through 
the inclusion of coupling modes between blade rows. It also enables the study of nonlinear effects of high 
amplitude blade vibration by the inclusion of higher harmonics of the fundamental blade vibration 
frequency. In the present work, the solver is applied to study in detail the aeroelastic characteristics of a 
transonic compressor stage. Various parameters were included in the study: number of coupling modes, 
blade row axial spacing, and operating speeds. Only the first vibration mode is considered with amplitude 
of oscillation in the linear range. Both aeroelastic stability (flutter) of rotor blade and unsteady loading 
on the stator are calculated. The study showed that for the stage considered, the rotor aerodynamic 
damping is not influenced by the presence of the stator even when the axial spacing is reduced by nearly 
25 percent. However, the study showed that blade row interaction effects become important for the 
unsteady loading on the stator when the axial spacing is reduced by the same amount. 

Introduction 

Turbomachinery designers try continually to improve the performance of aircraft engine components. 
However, the resulting complex blade shapes and high blade loadings, sometimes lead to undesirable 
unsteady aerodynamic and aeroelastic characteristics that affect the operability and durability of the 
engine components. In addition, the push towards reducing engine weight, leads to turbomachinery 
designs with fewer stages that are more closely spaced and have higher loading, which results in an 
increase in unsteady aerodynamic interactions between blade rows. Also, to reduce weight, integrally 
bladed disks (blisks) are increasingly used instead of the conventional inserted blades, reducing the 
amount of structural damping available to suppress the aeroelastic instabilities and forced response 
vibrations. 

Aeroelastic phenomena, which arise from the interaction between aerodynamic and structural 
dynamics, lead to blade vibrations. Blade vibrations can cause sudden blade failure as in flutter, or failure 
due to high cycle fatigue as in forced response. Flutter encountered unexpectedly during engine 
development and testing can result in long and costly delays in the development schedule due to the need 
to re-design and re-test the affected components and the engine. Forced response encountered by engines 
in service can cause High Cycle Fatigue (HCF) failures that lead to increased maintenance costs due to 
the increased need for inspection and repair. In-service engine failures can, and have, resulted in the loss 
of life. Also, such accidents are often followed by the grounding of an aircraft fleet, resulting in the 
disruption of scheduled service, increased operating costs, and ultimately in the loss of revenue. Therefore 
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it is highly important that the turbomachinery is designed to avoid any aeroelastic instabilities and 
vibrations. 

Current aeroelastic analysis of turbomachinery is focused on using time accurate computational fluid 
dynamic (CFD) approaches, i.e., high fidelity numerical techniques, for example, Reference 1. The 
methods and tools have resulted in predictions that showed good agreement with test data (Ref. 2). 
However the computational times required to run such codes are very high. Subsequently, research has 
focused on developing CFD based fast aeroelastic analysis models utilizing linearization. In these 
methods (Ref. 3) the governing unsteady nonlinear aerodynamic equations are written as steady non- 
linear equations and unsteady linearized equations. The unsteady linearized equations are obtained by 
linearizing the full non-linear unsteady equations about the nonlinear steady solution. These time- 
linearized formulations are fast and include steady loading effects. These methods have been applied to 
analyze turbomachinery flows (Ref. 4). However, this formulation neglects the dynamic nonlinear effects. 

To overcome the drawbacks in the methods mentioned above, harmonic balance techniques have 
been developed that can include dynamic nonlinear effects and at the same time can reduce time for 
solution (Refs. 5 and 6). In this approach, the unsteady flow is assumed to be temporally and spatially 
periodic, as in a turbomachinery flow. The flow variables are represented by a Fourier series in time with 
frequencies that are integer multiples of the original excitation frequency. The dependent variables are the 
coefficients of the Fourier series for each conservation variables. The Fourier series are then inserted into 
the unsteady fluid dynamics equations, and the resulting expressions are balanced by collecting terms 
with same frequency. For the equations to be satisfied, each frequency component must vanish 
independently. This results in a set of coupled complex partial differential equations, one for each 
frequency retained in the model. The method, which was later improved by storing the dependent 
variables at a number of sub-time levels over one period, has been applied to many configurations 
(Refs. 5 to 11). 

Most of the applications of the models mentioned above neglected the influence of neighboring blade 
rows. Analytical formulations have considered rigid blades to model the interaction effects and neglected 
blade vibration effects. An idea to include the neighboring blade rows including blade vibration was 
developed in Reference 12. This formulation used aerodynamic force coefficients obtained using linear 
unsteady cascade aerodynamic theory. Hall et al. (Refs. 13 and 14) extended this methodology to include 
CFD methods based on time-linearized formulations. Finally, the harmonic balance technique was 
combined with the multistage formulation, to include both multistage effects and nonlinear aerodynamic 
effects in the unsteady aerodynamic analysis (Refs. 15 to 18). 

In the present paper, the combined harmonic balance and multistage method, developed in 
References 13 to 18, is used to analyze the flutter and forced response characteristics of a compressor 
stage. The objectives of this paper are: (1) to apply the harmonic balance method with blade row 
interactions to obtain flutter and forced response characteristics of a transonic compressor stage for 
various operating conditions, (2) to study the effect of number of spinning modes on aerodynamic 
damping and unsteady loading, and (3) to investigate the effect of spacing of the blade rows on flutter and 
response. This is an assessment study for future application of the harmonic balance code to other 
multistage configurations that are being designed and tested at NASA. 

Formulation 

An aeroelastic model requires a consistent combination of a structural model with an unsteady 
aerodynamics model. The structural model used in the present analysis is obtained from a standard finite 
element analysis. The governing aeroelastic equations are recast in modal form. A brief description of the 
aeroelastic and the aerodynamic model is given below. See Reference 4 for the aeroelastic model 
formulation details. 
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Aeroelastic Model 


The aeroelastic equations of motion of s' h blade in any blade row can be written as 

= (i) 

where [A7 V ] and [A v ] are generalized mass and stiffness matrices, which are diagonal, { q s } is the 
generalized displacement vector, and [ A s ] is the motion-dependent generalized aerodynamic load matrix, 
and {AD S } is the motion-independent generalized aerodynamic load vector. The motion-dependent forces 
can cause flutter, and motion-independent forces result in forced response vibrations. The matrices \M S \, 
[A',.], and [A,] are of NMxNM size; {q s } and {AD S } are of NMx 1 size. 

For a tuned cascade (or rotor), in which all the blades are identical, the aeroelastic modes consist of 
individual blades vibrating with equal amplitudes with a fixed interblade phase angle between adjacent 
blades. Hence, the motion of the s h blade can be written as 

{q s }={qos}e W)t ={q a r}e ml e i °^ ( 2 ) 

Substituting in Equation (1), we obtain 

- oo 2 [M s ] {q ar + [K s ]{q ar Jgffa+M = [A r ]{q ar + {AD r (3) 

Since the blades are identical, the same equation is obtained for each blade and Equation (3) can be 
solved for N different values of the interblade phase angle given by 

(5 r 2 71// blades ? ^ ■ ~ ^blades /2 f / T Nblades /2 (4) 

where N b!ades is the number of blades in the cascade, and r is the nodal diameter (ND) index. 

Dropping the subscript s, since each blade is identical, and canceling out the exponential terms, 
Equation (3) can be written as 


- c« 2 [m}{y] + [A] {t} - [A r ]{7} = {AD,. } (5) 

Stability Calculation 

An eigenvalue approach is followed to investigate flutter stability. The flutter frequencies and flutter 
modes are obtained from the solution. For a stability calculation (flutter), the motion-independent forces 
{AD,) are set to zero, and Equation (5) can be written as 

- 0 2 [M]{7} + [A]{7}- [A r ]{Y}= {0} (6) 

Dividing Equation (6) with an assumed frequency, oo 2 

- (co/ C 05 )[at]{t} + [ [a] - [A,. ] ]{t}/(o^ = {0} (7) 

Rearranging, the equations can be written in the standard eigenvalue problem as: 

[fc]-y[e,-]]fr-Mo} (8) 

where 

fcHM - k ]]/® 2 ( 9 ) 
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and 


\Qr]=[M] 


( 10 ) 


y = (co/co 0 ) 2 (11) 

The solution of the above eigenvalue problem results in NM complex eigenvalues of the form 

i — -i yfy- p + i v (12) 

C0 o 

The real part of the eigenvalue ( p ) represents the damping ratio, and the imaginary part ( v ) 
represents the damped frequency; flutter occurs if p > 0 for any of the eigenvalues. For the tuned 
cascade, the stability of each phase angle mode is examined separately. The critical phase angle is 
identified as the one, which results in the lowest flutter speed. 

Response Calculation 

The aeroelastic response of the blades induced by wakes is calculated from Equation (5) as 

{Yr} = [[Pr]-l[Qr]V{AD r } (13) 

The amplitude of each blade is obtained by summing contributions of (T) from all phase angle modes. 

Aerodynamic Model 

The unsteady flow equations can be solved to obtain the elements of \A r ] and {AD,) using either a 
time domain approach or a frequency domain approach. Methods based on time domain approach are 
computationally intensive. Frequency domain approaches take out the terms associated with time by 
assuming a frequency for the solution, thereby reducing the computational time. In the current paper, the 
unsteady flow equations are solved using a frequency domain approach. Current frequency domain 
approaches are based on small amplitude oscillations, and assume either uniform steady flow (linear 
models), or a non-linear steady flow (time-linearized models). However, when the amplitude of 
oscillations is high, these models are not valid. A harmonic balance method can be developed that can 
still be in frequency domain and also can be used for high amplitude oscillations (Refs. 5 to 8). In the 
present formulation the harmonic balance method developed in Reference 5 is used. 

Harmonic Balance Method (Non-Linear Frequency Domain Model) 

In this approach, the conservation variables in the nonlinear unsteady governing equations are 
represented as a Fourier series in time with spatially varying coefficients, for example, 

p(x, y,t)=^R„ (x, y)eJ (ont 

n 

p u(x,y,t) = 'YjU,, (x, y)eJ wnt 

n 

pv(x, y, t) = 'Y J V n (x,y)ej ant 

n 
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pMx, y, 0 = 2] W n (x, y)eJ wnt 

n 

pe(x, y,t) = 'Y J E n (x, y)eJ ant ( 1 4) 


where, in principle, the summations are taken over all integer values of n. In practice, these series are 
truncated to a finite number of terms, —N<n< +N . Similar equations can be written for other 
conservation variables, and can be substituted into the Euler/Navier-Stokes equations. After substitution, 
the terms are grouped by frequency and require each frequency component to satisfy the conservation 
equations individually. Fourier components for non-negative n need to be stored. If A Harmonics are 
retained in the Fourier series representation of the flow, then 2/V+ 1 coefficients are stored for each flow 
variable (one for the zeroth harmonic or mean flow and 2 N for the real and imaginary parts of each 
remaining harmonic). This direct approach of substituting harmonics in governing equations has been 
shown to increase the computational time, and not be readily applicable to viscous flows because 
turbulence models tend to be quite complex and cannot be readily expressed in simple algebraic forms. 

To alleviate these problems, solutions at several sub-time levels that span a single time period are stored 
(Ref. 18). These sub-time level solutions are related to each other through the time derivative term, which 
is approximated by a pseudo-spectral operator in the governing equations, and through boundary 
conditions. The complex periodicity conditions connect the sub-time levels within a blade passage, and 
inter-row boundary conditions connect the solutions among the blade rows. The final discretized 
equations are mathematically steady, and the flowfield can be solved very efficiently using multi-grid 
acceleration and local time stepping techniques. The computational domain is reduced to a single blade 
passage, with complex periodic conditions applied along the periodic boundaries. The computational time 
then scales as the product of the number of sub-time levels and the time required for a single steady flow 
calculation. 

Coupling of Blade Rows 

In multistage turbomachinery, waves are created because of the unsteady loading on the blades. The 
waves from one blade row propagate upstream and downstream, and interact with the neighboring blade 
rows. The neighboring blade rows emit additional waves that propagate again upstream and downstream 
and interact with the blade rows adjacent to them. The waves are pressure waves, vorticity and entropy 
waves. This type of coupling is termed as aerodynamic coupling. If blades of one of the blade rows are 
vibrating, they emit additional waves that depend on their vibrations. The coupling due to blade vibrations 
is called dynamic coupling. In the present formulation, both the couplings are included by defining 
spinning modes, which are pressure/ entropy/ vorticity waves of different wave numbers (Refs. 12 to 18). 
A brief description is given here and additional details can be found in References 13 to 18. 

Consider a stage having a stator with B\ blades and a rotor with B 2 blades. Fet the original disturbance 
have k 0 nodal diameters (same as r in Eq. (4)) and frequency co () in the rotating frame reference. The 
acoustic, vortical and entropy waves from the stator blade row will interact with the rotor blade row, with 
the nodal diameter given by 


N = ko+ 1* B\+m* Bi 

with / and m taking on all integer values and frequencies given by 

oo' = coq + ffrA'o + m * Bj) for the stator 


(15a) 


(15b) 


and 
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oo" = coq - iT* l* B\ for the rotor 


(15c) 


Note that the frequencies in the stator frame of reference depend on the spinning mode index m associated 
with scattering of waves by the rotor, and the frequencies in the rotor frame of reference depend on the 
spinning mode index / associated with scattering of waves by the stator. At the beginning of the 
computation, the user specifies the spinning mode numbers, / and m to be used in the calculation. As an 
example, for one spinning mode, (l,m) = (0,0); for four spinning modes (l,m) = (0,0), (0,1), (1,0), (1,1); 
for nine modes (/,»?) = ( - 1 < / < 1,-1 < m < 1 ), and for 25 modes (/, w ) = ( - 2 < / < 2,-2 < m < 2 ), and so 
on. These mode numbers, together with the values of initial frequency, initial interblade phase angle, 
blade counts, and rotational speed determine the frequencies and the interblade phase angle of individual 
spinning modes. Unsteady solutions are generated for each of the spinning modes and combined to get the 
final unsteady solution. 

Theoretically, / and m can take on an infinite number of values, which generates an infinite number of 
spinning modes. As the number of spinning modes in the model is increased, the unsteady solution 
converges to a fixed solution or at least the solution of the fundamental mode converges. However, good 
estimates of the unsteady flowfield can be obtained using a relatively small number of spinning modes. The 
present solver calculates a number of intermediate flow field solutions for each blade row depending on the 
number of spinning modes. However, it is important to note that the number of intermediate solutions 
required to obtain the final unsteady flow field is far less than the number of spinning modes selected, since 
many spinning modes have same frequencies and interblade phase angles differing only by 2n . 

Finally, the approach can be extended to include higher harmonics of vibration amplitude as 

N = n* ko + 1 * B\ + m* Bi for nodal diameters (16a) 

oo' = n * coq + Q(k o + m * B 2 ) frequency for the stator ( 1 6b) 

co" = «*cc>o - Ci*l*Bi frequency for the rotor (16c) 

where n can be 1,2,.. .NH, harmonics of vibration. Note that selected values of ( n,l,m ) are given as the 
input. Also, note that n= 1 is implied in Equations (15). 

Results 

The configuration selected for application is designated as R67-S67 compressor stage. Experimental 
data is available in References 19 and 20. Rotor 67 is a transonic axial flow fan rotor tested by NASA in a 
rotor alone mode (Ref. 19). This rotor was designed with 22 blades with a design relative inlet Mach 
number at the rotor tip of 1.38. The design rotational speed is 16,043 rpm, and the total pressure ratio is 
1.63 at a mass flow of 33.25 kg/sec. To extend to a single stage, a blade row of stators (Stator 67) with 34 
blades was placed downstream of the rotor row (Ref. 20). The axial distance from the trailing edge of the 
rotor to the leading edge of the stator is 85 percent of rotor axial chord at mid-span. This configuration 
was studied in Reference 1 6 for one operating point. In the present work, the study is extended to multiple 
parameters, including spinning modes, operating speeds, and aerodynamic damping calculations. 

Calculations were carried out at the design rotational speed of 100 percent (16,043 lpm) and part 
speed of 75 percent (12,032 rpm). It was assumed that changes in blade geometry due to changes in 
rotational speed were not significant and therefore a single geometry (at 1 00 percent speed) was used for 
computations at all speeds and operating conditions. The computational grid used, generated using 
commercial software, is shown in Figure 1; the grid size is 193x25x33 for the block that wraps around the 
blade airfoil (O-grid) with 1 93 grid points around the airfoil; 25 grid points in the circumferential 
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(a) Meridional view 


(b) Azimuthal view 


Figure 1. — Rotor-stator configuration and computational grids. 

direction, and 33 grid points in the spanwise direction. The grid blocks (H-grid) in the inlet and exit 
sections are each 17x25x33 with 17 points in the streamwise direction, 25 grid points in the 
circumferential direction, and 33 grid points in the spanwise direction. The tip clearance was modeled 
with a grid with 5 points between the blade tip and casing. The same type of grids was used for the stator. 
However, there was no clearance grid for the stator. Figure 1 shows the rotor-stator arrangement. 

Uniform inlet flow conditions (constant total pressure and temperature) were used in the 
computations. The exit flow conditions consisted of a static pressure prescribed at the casing and 
calculated at other radial locations using radial equilibrium. All solutions were obtained using a single 
blade passage for each blade row. 

Steady Solution at 100 percent Speed 

The first step in the validation of the code is to simulate the operating map of the compressor stage. 
Steady calculations were performed for various values of the imposed non-dimensional back-pressure, 
p exit, at the exit boundary for 100 percent rotational speed. Note that for steady calculation, all the 
spinning mode parameters, Equation (16) are made equal to 0; i.e. ( n,l,m ) = (0,0,0). Figure 2 shows 
the convergence of the steady computations as a plot of non-dimensional torque with iteration counter 
(100 iterations per counter) for the rotor. As can be noted from Figure 2, excellent convergence was 
obtained for all operating conditions. Similar convergence was obtained for the stator and for other 
computations at 75 percent speed presented later. All steady results presented here are from well- 
converged solutions. 

Figure 3 shows the computed operating map with 100 percent speed line. Solutions were obtained for 
various values of p exit to obtain the complete operating characteristic (speed line) for pressure ratios that 
include both stall and choke sides. Note that Reference 16 previously presented steady flow at a single 
operating point on the speed line to validate the code. In the present work, the flow solutions were post- 
processed by a routine that calculates pressure ratio and mass flow rate. 

Figure 3(a) shows the total pressure ratio versus mass flow rate calculated, between the inlet and the 
stage exit, by the in-house code PERF, and is designated as plp3; published experimental measurements 
from Reference 20 are also included for comparison. It can be seen that the present computational results 
are shifted towards a lower mass flow rate (approximately 1 kg/s) and a lower pressure ratio 
(approximately 0.1) demonstrating that present code predicts the steady flow reasonably well. 
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The predicted adiabatic efficiency is shown in Figure 4. Digitized test values from Reference 20 are 
also depicted in the figure. It is noted that the calculated efficiency is lower than the measured value for 
all the mass flow rates. The calculated peak efficiency occurs at a mass flow rate of 32.75 kg/sec 
compared to the test value of 34.0 kg/sec, a difference of 3.7 percent in the mass flow rate. However, the 
sudden jump in the efficiency at 32.75 kg/sec requires further study. 



0 10 20 30 40 50 

Iteration Counter 

Figure 2. — Convergence of the steady solution for different values of back- 
pressure, p_exit. 
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Figure 3(a). — Stage operating map. 
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Figure 3(b). — Stage and rotor operating map. 
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Figure 4. — Stage adiabatic efficiency. 


Unsteady Calculations: Flutter 

Unsteady computations were performed using the harmonic balance code to predict the aeroelastic 
stability of the rotor. The first vibrating rotor mode is considered in this study. As shown in Figure 5, the 
mode is a predominantly bending mode with a twisting component at a frequency of 554 Flz. The 
prescribed blade vibration amplitude is in the linear range. All solutions were obtained using a single 
blade passage for each blade row. 
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100 Percent Speed 
Fundamental Spinning Mode 

Results of unsteady flow calculations, obtained using the fundamental spinning mode (, n,l,m ) = (1,0,0) 
as per Eq. (16)), are presented in this section. These computations were for an operating point for a non- 
dimensional back-pressure p exit of 1.0. The unsteady solution was started from the steady condition. 
Several nodal diameters ( ND ) covering the entire possible range, Equation (4), were considered. For the 
22 blade rotor considered here, the nodal diameters vary from -1 1 to 1 1 as per Equation (4). Figure 6 
shows the convergence of the unsteady computations as a plot of non-dimensional generalized force 
(imaginary part) with iteration counter (100 iterations per counter) for various nodal diameters. 

Figure 7 shows the unsteady generalized forces on the rotor calculated for various nodal diameters for 
the same back-pressure, p exit, of 1.0. Real values are related to the frequency, and imaginary values are 
related to the damping and stability of the system. The computed values compare well with the published 
results. For a single vibration mode, the stability is determined by the imaginary part of the generalized 
force, and Figure 7 shows that the stability varies significantly with nodal diameter, it can be seen that the 
imaginary part of the generalized force has the lowest absolute value for ND=0, 1, and 2. 

The damping of the system is obtained from the generalized forces by solving the eigenvalue problem 
given in Equation (8). The damping obtained from Equation (12) for various nodal diameters is plotted in 
Figure 8. Note that the damping values presented in the figures are opposite in sign to the damping ratio 
(real part of the eigenvalue, p ) in Equation (12). Structural damping is not considered in the calculations; 
therefore, the calculated damping from Equation (12) is entirely aerodynamic damping. It can be seen that 
the damping is the lowest for ND=0, 1, and 2, and varies from 0.5 to 3.5 percent for the range of nodal 
diameters. 



Figure 5. — First mode shape of 
the rotor blade (554 Flz). 
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Figure 6. — Convergence of the unsteady generalized force (imaginary 
part). 


NASA/TM— 2010-216240 


10 


20 



-11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 

Nodal Diameter 

Figure 7. — Rotor generalized force (real and imaginary parts); p_exit=1.0. 
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Figure 8. — Rotor aerodynamic damping at 100 percent speed; p_exit=1.0. 


Additional Spinning Modes 

In addition to the fundamental spinning mode, ( n,l,m ) = (0,0,0), additional spinning modes, are 
included in the analysis to study their effect on improving the accuracy of the unsteady flow solution. The 
analysis is still done at small amplitude oscillations by including only one harmonic («= 1). The modes 
selected are ( n,l,m ) = (1,0,0), (1,0,1), (1,1,0) and (1,1,1) as per Equations (16). Figure 9 shows the 
generalized forces calculated with the fundamental mode (designated nmod2) and extra spinning modes 
(designated nmod=5). It can be seen that the addition of more spinning modes in the solution did not 
affect the generalized force. Additional solutions were also obtained by including more modes (nine 
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modes), but they did not show any appreciable change in the generalized force. These results are therefore 
not included here. The results in Figure 9 clearly show that for this configuration, at the conditions 
considered, the stator is not influencing the rotor in any significant way. Therefore, in the following 
sections, only the results obtained using two spinning modes are presented. 

Effect of Back-Pressure, pexit 

Figure 10 shows the aerodynamic damping for back-pressure of 1.0 and 1.055. It can be noted that the 
aerodynamic damping calculated for all nodal diameters is positive. The stability varies slightly with 
back-pressure. Further, with the change in back-pressure from peak-efficiency towards stall, it can be 
seen that the damping decreases and the graph moves towards the zero damping line. The nodal diameters 
corresponding to lowest damping remain unchanged (ND= 0, 1, 2). However, the reduction in damping 
with back-pressure is marginal. For the last computed point near stall (p_exit=1.055), the damping 
remains positive, indicating that no flutter is expected for the 1 00 percent speed. This result is consistent 
with the rig test, in which no flutter was observed. 
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Figure 9. — Effect of spinning modes on rotor generalized forces. 
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Figure 10. — Effect of back-pressure on damping. 



Blade Row Spacing Effects 

The results presented in the preceding sections showed that the blade row interaction effects for this 
configuration are marginal. Hence, to increase the interaction between blade rows, the axial spacing 
between the rotor and stator is reduced from 85 percent of rotor chord to 65 percent of rotor chord 
(approximately 25 percent reduction in axial spacing). The revised operating map for this configuration is 
shown in Figure 1 1 with the results for the original spacing (85 percent rotor chord) designated as spl, 
and the results for the reduced spacing (65 percent rotor chord) designated as sp2. It can be seen that the 
rotor speed line moved up towards higher pressure ratios for the same mass flow rate. This result is 
consistent with the expected effect of a reduction in blade row axial spacing. 

Figure 12 shows the generalized force obtained when the spacing between the rotor and stator is 
reduced by 25 percent. Only the fundamental mode is included in these calculations since previous results 
indicated no significant effect of the additional spinning modes. In Figure 12, it can be seen that the 
generalized forces varied somewhat for some nodal diameters such as ND=S and 1 1 . However, for the 
nodal diameters with the least damping, there is minimal change in the unsteady generalized force on the 
rotor even with the 25 percent reduction in spacing. It is seen that the stator has minimal influence on the 
rotor unsteady loading. It is expected that the stator may be more influenced by the rotor than vice versa 
and this possibility will be examined later by comparing the unsteady loading on the stator for both the 
spacings. 


75 Percent Speed 

To examine the possibility of part-speed flutter, computations were carried out at 75 percent 
rotational speed (12,032 rpm). Figures 13 and 14 show the operating map with 75 and 100 percent speed 
lines. The 75 percent speed line is shifted to lower mass flow rate and pressure ratio, as expected. The 
drop in the pressure ratio is about 0.4, similar to published data for similar configurations. The reduction 
in mass flow rate is about 8 kg/sec, which is also similar to published data for similar configurations. The 
peak efficiency occurs at a pressure ratio of 25.65 kg/sec for a pressure ratio of 1.290. 
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Figure 13. — Rotor operating map. 


Figure 15 shows the variation of aerodynamic damping with nodal diameter at 75 percent rotational 
speed and for the original blade row spacing (spl). The vibration frequency is 479 Flz for this speed. The 
figure shows the calculated aerodynamic damping for two values of back-pressure (pexit). The lower 
value (0.85) is near peak efficiency, and the higher value (0.92) is near stall. It can be seen that the 
damping of the rotor is reduced at 75 percent speed as compared to the damping at 100 percent speed 
(Fig. 10), but only for high nodal diameters. At the low nodal diameters, the damping levels are nearly 
unchanged between 100 and 75 percent speeds. 

To examine the stability characteristics at near-stall conditions, the aerodynamic damping at 
100 percent speed and 75 percent speed is compared in Figure 16. The corresponding operating points 
near stall are p_exit=T.055 for 100 percent speed and p_exit=0.92 for 75 percent speed. It can be seen that 
the lowest value of damping near nodal diameters (0,1, and 2) remained the same, even though the highest 
value of the damping dropped from 3.5 percent at 100 percent speed to 2.5 percent at 75 percent speed 
line. Thus, for this configuration, the lowest aerodynamic damping levels are nearly the same at both the 
design and part speeds. Further, the nodal diameters corresponding to the lowest damping are also the 
same for design and part speeds. 
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Figure 14. — Rotor Efficiency map. 
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Figure 15. — Effect of back-pressure (p_exit) on aerodynamic damping at 
75 percent speed. 
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Figure 16. — Effect of speed on aerodynamic damping for back-pressures 
near stall. 


Unsteady Calculations: Response 

For the forced response vibrations of the stator blades, calculations were carried without blade 
vibrations. The focus is on the forcing function acting on the stator blades due to the presence of the rotor. 
A single blade passage is used for each blade row in this response study similar to the flutter calculations 
described earlier. For most of the calculations, four spinning modes are selected for the response solution: 
(: n,l,m ) = (0,0,0), (0,1,1), (0,1,0) and (0,1,1) as per Equation (16). This selection results in single blade 
passing frequency. Note that n = 0 for all the spinning modes since blade vibration is not included in 
response calculations. Additional calculations were also performed with 9 and 16 modes, which resulted 
in two and three discrete non-zero frequencies, following Reference 2 1 . The frequencies of the forcing 
function will be the blade passing frequency and its higher harmonics. For this case, the frequency will be 
5,882 Hz, which is about 10 times the vibration frequency used in the rotor flutter study. 

Figure 17 shows the mean (time-averaged) pressure distribution on the blade surface at 90 percent 
span of the stator blade row obtained with 4 and 9 modes, and for the original blade row spacing (spl). 
Also shown in the figure is the pressure distribution for the multiblade row steady solution 
(n,l,m)=( 0,0,0). The mean surface pressure obtained from the unsteady computations is slightly different 
than the one obtained from the steady computation, although the differences are small. The results from 
solutions with 4 and 9 modes are nearly identical, indicating that the mean pressure solution has 
converged with respect to the number of modes. 

Figure 18 shows the instantaneous unsteady pressure contours at 90 percent span from 4 mode 
solution. Although the actual computations are performed over a single blade passage in each blade row, 
multiple passages are included in Figure 18 for better visualization of the flow features. The contour plots 
show no significant flow features that may cause a strong interaction between the blade rows. 

Figure 19 shows the real and imaginary parts of the first harmonic of the unsteady pressure at 
90 percent span of the stator. Small differences are seen between the solutions for 4 and 9 modes. The 
solution with 9 modes shows slight improvement in the spatial resolution of the pressure distribution, as 
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can be seen in Figure 19(a) from leading edge to 40 percent chord. Based on these representative results, 
it is inferred that four modes are enough to capture the primary unsteady effect for the present study. 
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Figure 17. — Mean surface pressure distribution at 90 percent span of the 
stator blade, p_exit=1.0. 
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Figure 18. — Instantaneous pressure contours at 90 percent span on 
stator blade; 4 mode solution, p_exit=1.0. 
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Figure 19. — Unsteady first harmonic surface pressure distribution at 
90 percent span of the stator blade, p_exit=1 .0 
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Figure 20 shows unsteady surface pressure distribution at 90 percent span of the stator blade row for 
the original blade row spacing (spl) and reduced blade row spacing (sp2). Four modes are used for both 
these solutions. There are notable differences in the unsteady loading for the two spacings considered. 
Based on these results, it may be inferred that the unsteady forcing function on the stator is changed by 
the reduction in the axial spacing between the two blade rows. A comparison of the generalized force for 
the two axial spacings will be required to quantify the impact of spacing on stator blade response. 



x/c, Chordwise location 

(a) real part 
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(b) imaginary part 

Figure 20. — Unsteady first harmonic surface pressure distribution at 
90 percent span of the stator blade for two axial spacings; p_exit=1 .0 
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Conclusion 


A harmonic balance method has been applied to analyze flutter and response of a compressor rotor- 
stator configuration. The unsteady aerodynamic analysis includes blade row interaction effects through 
the inclusion of coupling modes between blade rows. The compressor configuration analyzed was 
previously tested at NASA but blade vibrations were not measured or studied in the experiment. The 
steady calculations were done for multiple operating points and rotational speeds than those published in 
literature. The present calculations predicted the steady loading close to that of the design and the value 
measured in experiment. 

Unsteady flutter calculations were performed for the rotor-stator configuration. In this study, the first 
vibration mode of the rotor with amplitude of vibrations in the linear range was considered. Aerodynamic 
damping results showed a very weak variation with back-pressure or mass flow rate. Detailed calculations 
indicated that the rotor was stable for the rotational speeds and mass flow rates considered. Increasing the 
number of coupling modes used in the unsteady solution did not change the results obtained using only 
the fundamental coupling mode. A reduction of the axial spacing between rotor and stator by 25 percent 
showed no significant effect on the aerodynamic damping results. For all the cases studied, the low nodal 
diameters of ND= 0, 1 , 2 were found to have the lowest aerodynamic damping. 

The unsteady loading on the stator was also calculated using the present harmonic balance method. 
Solutions with 4, 9, and 16 coupling modes were obtained and compared. It was noted that the mean 
loading from the unsteady solution was slightly different from loading from the solution with no coupling 
modes. Small differences were seen between the solutions for 4 and 9 modes. The solution with 9 modes 
shows slight improvement in the spatial resolution of the pressure distribution. Based on these 
representative results, it is inferred that four modes are enough to capture the primary unsteady effect for 
the present study. Reduction of the axial spacing between rotor and stator by 25 percent showed 
considerable effect on the unsteady loading on the stator. A comparison of the generalized force for the 
two axial spacings is planned to quantify the impact of spacing on stator blade response. 

Follow-on work is planned to include other modes of rotor vibration for the R67-S67 compressor 
configuration and response calculations of the stator blades. Further application and validation of the 
harmonic balance method is planned for the unsteady analysis of multistage turbomachinery. 
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